In [1]:

    

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import astropy.units as u
from astropy import time


    

In [2]:

    

from poliastro import iod
from poliastro.plotting import plot
from poliastro.bodies import Sun, Earth
from poliastro.twobody import State
from poliastro import ephem

from jplephem.spk import SPK
ephem.download_kernel("de421")


    

    




File de421.bsp already exists under /home/juanlu/.poliastro

In [3]:

    

r = [-6045, -3490, 2500] * u.km
v = [-3.457, 6.618, 2.533] * u.km / u.s

ss = State.from_vectors(Earth, r, v)

with plt.style.context('pybonacci'):
    plot(ss)


    

In [4]:

    

epoch = time.Time("2015-06-21 16:35")

r_, v_ = ephem.planet_ephem(ephem.EARTH, epoch)


    

In [5]:

    

r_


    

    
Out[5]:





$[-61689.755,~-1.394555 \times 10^{8},~-60480590] \; \mathrm{km}$

In [6]:

    

v_.to(u.km / u.s)


    

    
Out[6]:





$[29.308324,~-0.19527688,~-0.084878762] \; \mathrm{\frac{km}{s}}$

In [7]:

    

date_launch = time.Time('2011-11-26 15:02', scale='utc')
date_arrival = time.Time('2012-08-06 05:17', scale='utc')
tof = date_arrival - date_launch
r0, _ = ephem.planet_ephem(ephem.EARTH, date_launch)
r, _ = ephem.planet_ephem(ephem.MARS, date_arrival)

(v0, v), = iod.lambert(Sun.k, r0, r, tof)


    

In [8]:

    

v0


    

    
Out[8]:





$[-29.29151,~14.533265,~5.4169134] \; \mathrm{\frac{km}{s}}$

In [9]:

    

v


    

    
Out[9]:





$[17.615499,~-10.998307,~-4.2079605] \; \mathrm{\frac{km}{s}}$

In [10]:

    

def go_to_mars(offset=500., tof_=6000.):
    # Initial data
    N = 50

    date_launch = time.Time('2016-03-14 09:31', scale='utc') + ((offset - 500.) * u.day)
    date_arrival = time.Time('2016-10-19 16:00', scale='utc') + ((offset - 500.) * u.day)
    tof = tof_ * u.h

    # Calculate vector of times from launch and arrival Julian days
    jd_launch = date_launch.jd
    jd_arrival = jd_launch + tof.to(u.day).value
    jd_vec = np.linspace(jd_launch, jd_arrival, num=N)

    times_vector = time.Time(jd_vec, format='jd')
    rr_earth, vv_earth = ephem.planet_ephem(ephem.EARTH, times_vector)
    rr_mars, vv_mars = ephem.planet_ephem(ephem.MARS, times_vector)
    # Compute the transfer orbit!
    r0 = rr_earth[:, 0]
    rf = rr_mars[:, -1]

    (va, vb), = iod.lambert(Sun.k, r0, rf, tof)

    ss0_trans = State.from_vectors(Sun, r0, va, date_launch)
    ssf_trans = State.from_vectors(Sun, rf, vb, date_arrival)
    # Extract whole orbit of Earth, Mars and transfer (for plotting)
    rr_trans = np.zeros_like(rr_earth)
    rr_trans[:, 0] = r0
    for ii in range(1, len(jd_vec)):
        tof = (jd_vec[ii] - jd_vec[0]) * u.day
        rr_trans[:, ii] = ss0_trans.propagate(tof).r

    # Better compute backwards
    jd_init = (date_arrival - 1 * u.year).jd
    jd_vec_rest = np.linspace(jd_init, jd_launch, num=N)

    times_rest = time.Time(jd_vec_rest, format='jd')
    rr_earth_rest, _ = ephem.planet_ephem(ephem.EARTH, times_rest)
    rr_mars_rest, _ = ephem.planet_ephem(ephem.MARS, times_rest)
    # Plot figure
    # To add arrows:
    # https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/streamplot.py#L140
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')

    def plot_body(ax, r, color, size, border=False, **kwargs):
        """Plots body in axes object.

        """
        return ax.plot(*r[:, None], marker='o', color=color, ms=size, mew=int(border), **kwargs)

    # I like color
    color_earth0 = '#3d4cd5'
    color_earthf = '#525fd5'
    color_mars0 = '#ec3941'
    color_marsf = '#ec1f28'
    color_sun = '#ffcc00'
    color_orbit = '#888888'
    color_trans = '#444444'

    # Plotting orbits is easy!
    ax.plot(*rr_earth.to(u.km).value, color=color_earth0)
    ax.plot(*rr_mars.to(u.km).value, color=color_mars0)
    ax.plot(*rr_trans.to(u.km).value, color=color_trans)

    ax.plot(*rr_earth_rest.to(u.km).value, ls='--', color=color_orbit)
    ax.plot(*rr_mars_rest.to(u.km).value, ls='--', color=color_orbit)

    # But plotting planets feels even magical!
    plot_body(ax, np.zeros(3), color_sun, 16)

    plot_body(ax, r0.to(u.km).value, color_earth0, 8)
    plot_body(ax, rr_earth[:, -1].to(u.km).value, color_earthf, 8)

    plot_body(ax, rr_mars[:, 0].to(u.km).value, color_mars0, 8)
    plot_body(ax, rf.to(u.km).value, color_marsf, 8)

    # Add some text
    ax.text(-0.75e8, -3.5e8, -1.5e8, "ExoMars mission:\nfrom Earth to Mars",
            size=20, ha='center', va='center', bbox={"pad": 30, "lw": 0, "fc": "w"})
    ax.text(r0[0].to(u.km).value * 2.4, r0[1].to(u.km).value * 0.4, r0[2].to(u.km).value * 1.25,
            "Earth at launch\n({})".format(date_launch.to_datetime().strftime("%d %b")),
            ha="left", va="bottom", backgroundcolor='#ffffff')
    ax.text(rf[0].to(u.km).value * 1.1, rf[1].to(u.km).value * 1.1, rf[2].to(u.km).value,
            "Mars at arrival\n({})".format(date_arrival.to_datetime().strftime("%d %b")),
            ha="left", va="top", backgroundcolor='#ffffff')
    ax.text(-1.9e8, 8e7, 1e8, "Transfer\norbit", ha="right", va="center", backgroundcolor='#ffffff')

    # Tune axes
    ax.set_xlim(-3e8, 3e8)
    ax.set_ylim(-3e8, 3e8)
    ax.set_zlim(-3e8, 3e8)

    # And finally!
    ax.view_init(30, 260)
    plt.show()
    #fig.savefig("trans_30_260.png", bbox_inches='tight')
    #return fig, ax


    

In [11]:

    

go_to_mars()


    

In [12]:

    

%matplotlib inline
from ipywidgets import interactive
from IPython.display import display

w = interactive(go_to_mars, offset=(0., 1000.), tof_=(100., 12000.))
display(w)